# 03_PSM_Forced_Choice.R
# Created: 2021-05-03 Taka-aki Asano
# Last Modified: 2024-02-16

# package
require("readr")
require("dplyr")
require("tidyr")
require("ggplot2")
require("ggpubr")
require("cjoint")
require("cregg")
theme_set(theme_minimal(base_size = 12))


# load dataset
data <- read_csv("HLWconjoint20210330.csv")


# processing of data
data2 <- select(data, 
                no, time, rl, group, age, psm01, psm02, psm03, psm04, 
                psm05, psm06, psm07, psm08, psm09, psm10, psm11, 
                psm12, psm13, psm14, psm15, psm16, hlw_su, hlw_fsu, 
                hlw_dis, hlw_pol, hlw_so, hlw_ad, hlw_pop, hlw_inc)
data2$hlw_fsu <- ifelse(data2$hlw_fsu == data2$rl, 1, 0)
data2 <- mutate(data2, 
                `Distance` = recode(hlw_dis, `1` = "less than 1 km", 
                                    `2` = "1-5 km", `3` = "5-20 km", 
                                    `4` = "20 km or more"), 
                `Politics` = recode(hlw_pol, `1` = "job creation", 
                                    `2` = "public services", 
                                    `3` = "local economy", 
                                    `4` = "society"), 
                `Social` = recode(hlw_so, `1` = "70% agreed, 30% disagreed", 
                                  `2` = "50% agreed, 50% disagreed", 
                                  `3` = "30% agreed, 70% disagreed"), 
                `Administration` = recode(hlw_ad, `1` = "none", 
                                          `2` = "childbirth and childcare", 
                                          `3` = "elderly care", `4` = "education", 
                                          `5` = "public facility", `6` = "water charge"), 
                `Population` = recode(hlw_pop, `1` = "19,000", 
                                      `2` = "20,000", `3` = "21,000"), 
                `Income` = recode(hlw_inc, `1` = "3.8 million yen", 
                                  `2` = "4 million yen", `3` = "4.2 million yen"))
data2$Distance <- factor(
  data2$Distance, 
  levels = c("less than 1 km", "1-5 km", "5-20 km", "20 km or more")
)
data2$Politics <- factor(
  data2$Politics, 
  levels = c("society", "job creation", "public services", "local economy")
)
data2$Social <- factor(
  data2$Social, 
  levels = c("30% agreed, 70% disagreed", "50% agreed, 50% disagreed", "70% agreed, 30% disagreed")
)
data2$Administration <- factor(
  data2$Administration, 
  levels = c("none", "childbirth and childcare", "elderly care", 
             "education", "public facility", "water charge")
)
data2$Population <- factor(
  data2$Population, 
  levels = c("19,000", "20,000", "21,000")
)
data2$Income <- factor(
  data2$Income, 
  levels = c("3.8 million yen", "4 million yen", "4.2 million yen")
)


# group
## PSM
data2$PSM <- rowSums(data2[,paste0("psm", c("01", "02", "03", "04", "05", "06", 
                                            "07", "08", "09", "10", "11", "12", 
                                            "13", "14", "15", "16"))])
summary(data2$PSM)
data2$PSM_Group <- ifelse(data2$PSM >= median(data2$PSM), "HPSM", "LPSM")
data2$PSM_Group <- factor(data2$PSM_Group, levels = c("LPSM", "HPSM"))
aggregate(PSM ~ PSM_Group, data2, mean)

## APS
data2$APS <- rowSums(data2[,paste0("psm", c("01", "03", "12", "15"))])
summary(data2$APS)
data2$APS_Group <- ifelse(data2$APS >= median(data2$APS), "HAPS", "LAPS")
data2$APS_Group <- factor(data2$APS_Group, levels = c("LAPS", "HAPS"))
aggregate(APS ~ APS_Group, data2, mean)

## CPV
data2$CPV <- rowSums(data2[,paste0("psm", c("04", "06", "07", "08"))])
summary(data2$CPV)
data2$CPV_Group <- ifelse(data2$CPV >= median(data2$CPV), "HCPV", "LCPV")
data2$CPV_Group <- factor(data2$CPV_Group, levels = c("LCPV", "HCPV"))
aggregate(CPV ~ CPV_Group, data2, mean)

## COM
data2$COM <- rowSums(data2[,paste0("psm", c("02", "09", "10", "16"))])
summary(data2$COM)
data2$COM_Group <- ifelse(data2$COM >= median(data2$COM), "HCOM", "LCOM")
data2$COM_Group <- factor(data2$COM_Group, levels = c("LCOM", "HCOM"))
aggregate(COM ~ COM_Group, data2, mean)

## SS
data2$SS <- rowSums(data2[,paste0("psm", c("05", "11", "13", "14"))])
summary(data2$SS)
data2$SS_Group <- ifelse(data2$SS >= median(data2$SS), "HSS", "LSS")
data2$SS_Group <- factor(data2$SS_Group, levels = c("LSS", "HSS"))
aggregate(SS ~ SS_Group, data2, mean)


# PSM
## AMCEs
res1 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ PSM_Group, level_order = "descending"
)
res1$level <- as.character(res1$level)
res1 <- rbind(
  res1, 
  c("HPSM", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HPSM")
)
res1$level <- factor(
  res1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res1$estimate <- as.numeric(res1$estimate)
res1$std.error <- as.numeric(res1$std.error)
res1$z <- as.numeric(res1$z)
res1$p <- as.numeric(res1$p)
res1$lower <- as.numeric(res1$lower)
res1$upper <- as.numeric(res1$upper)
res1.plot <- ggplot(
  res1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = PSM_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res1.plot)

diff1 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ PSM_Group, level_order = "descending"
)
diff1$level <- as.character(diff1$level)
diff1 <- rbind(
  diff1, 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Politics", "(Political attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Social", "(Social attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Population", "(Population attribute)", "HPSM", NA, NA, NA, NA, NA, NA), 
  c("HPSM - LPSM", "hlw_fsu", "amce", "Income", "(Average income attribute)", "HPSM", NA, NA, NA, NA, NA, NA)
)
diff1$level <- factor(
  diff1$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff1$estimate <- as.numeric(diff1$estimate)
diff1$std.error <- as.numeric(diff1$std.error)
diff1$z <- as.numeric(diff1$z)
diff1$p <- as.numeric(diff1$p)
diff1$lower <- as.numeric(diff1$lower)
diff1$upper <- as.numeric(diff1$upper)
diff1.plot <- ggplot(
  diff1, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff1.plot)

## MMs
res2 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ PSM_Group, level_order = "descending"
)
res2$level <- as.character(res2$level)
res2 <- rbind(
  res2, 
  c("HPSM", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HPSM")
)
res2$level <- factor(
  res2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res2$estimate <- as.numeric(res2$estimate)
res2$std.error <- as.numeric(res2$std.error)
res2$z <- as.numeric(res2$z)
res2$p <- as.numeric(res2$p)
res2$lower <- as.numeric(res2$lower)
res2$upper <- as.numeric(res2$upper)
res2.plot <- ggplot(
  res2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = PSM_Group)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res2.plot)

diff2 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ PSM_Group, level_order = "descending"
)
diff2$level <- as.character(diff2$level)
diff2 <- rbind(
  diff2, 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HPSM"), 
  c("HPSM - LPSM", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HPSM")
)
diff2$level <- factor(
  diff2$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff2$estimate <- as.numeric(diff2$estimate)
diff2$std.error <- as.numeric(diff2$std.error)
diff2$z <- as.numeric(diff2$z)
diff2$p <- as.numeric(diff2$p)
diff2$lower <- as.numeric(diff2$lower)
diff2$upper <- as.numeric(diff2$upper)
diff2.plot <- ggplot(
  diff2, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff2.plot)


# APS
## AMCEs
res3 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ APS_Group, level_order = "descending"
)
res3$level <- as.character(res3$level)
res3 <- rbind(
  res3, 
  c("HAPS", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HAPS")
)
res3$level <- factor(
  res3$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res3$estimate <- as.numeric(res3$estimate)
res3$std.error <- as.numeric(res3$std.error)
res3$z <- as.numeric(res3$z)
res3$p <- as.numeric(res3$p)
res3$lower <- as.numeric(res3$lower)
res3$upper <- as.numeric(res3$upper)
res3.plot <- ggplot(
  res3, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = APS_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res3.plot)

diff3 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ APS_Group, level_order = "descending"
)
diff3$level <- as.character(diff3$level)
diff3 <- rbind(
  diff3, 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Politics", "(Political attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Social", "(Social attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Population", "(Population attribute)", "HAPS", NA, NA, NA, NA, NA, NA), 
  c("HAPS - LAPS", "hlw_fsu", "amce", "Income", "(Average income attribute)", "HAPS", NA, NA, NA, NA, NA, NA)
)
diff3$level <- factor(
  diff3$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff3$estimate <- as.numeric(diff3$estimate)
diff3$std.error <- as.numeric(diff3$std.error)
diff3$z <- as.numeric(diff3$z)
diff3$p <- as.numeric(diff3$p)
diff3$lower <- as.numeric(diff3$lower)
diff3$upper <- as.numeric(diff3$upper)
diff3.plot <- ggplot(
  diff3, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff3.plot)

## MMs
res4 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ APS_Group, level_order = "descending"
)
res4$level <- as.character(res4$level)
res4 <- rbind(
  res4, 
  c("HAPS", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HAPS")
)
res4$level <- factor(
  res4$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res4$estimate <- as.numeric(res4$estimate)
res4$std.error <- as.numeric(res4$std.error)
res4$z <- as.numeric(res4$z)
res4$p <- as.numeric(res4$p)
res4$lower <- as.numeric(res4$lower)
res4$upper <- as.numeric(res4$upper)
res4.plot <- ggplot(
  res4, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = APS_Group)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res4.plot)

diff4 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ APS_Group, level_order = "descending"
)
diff4$level <- as.character(diff4$level)
diff4 <- rbind(
  diff4, 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HAPS"), 
  c("HAPS - LAPS", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HAPS")
)
diff4$level <- factor(
  diff4$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff4$estimate <- as.numeric(diff4$estimate)
diff4$std.error <- as.numeric(diff4$std.error)
diff4$z <- as.numeric(diff4$z)
diff4$p <- as.numeric(diff4$p)
diff4$lower <- as.numeric(diff4$lower)
diff4$upper <- as.numeric(diff4$upper)
diff4.plot <- ggplot(
  diff4, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff4.plot)


# CPV
## AMCEs
res5 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ CPV_Group, level_order = "descending"
)
res5$level <- as.character(res5$level)
res5 <- rbind(
  res5, 
  c("HCPV", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCPV")
)
res5$level <- factor(
  res5$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res5$estimate <- as.numeric(res5$estimate)
res5$std.error <- as.numeric(res5$std.error)
res5$z <- as.numeric(res5$z)
res5$p <- as.numeric(res5$p)
res5$lower <- as.numeric(res5$lower)
res5$upper <- as.numeric(res5$upper)
res5.plot <- ggplot(
  res5, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = CPV_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res5.plot)

diff5 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ CPV_Group, level_order = "descending"
)
diff5$level <- as.character(diff5$level)
diff5 <- rbind(
  diff5, 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Politics", "(Political attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Social", "(Social attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Population", "(Population attribute)", "HCPV", NA, NA, NA, NA, NA, NA), 
  c("HCPV - LCPV", "hlw_fsu", "amce", "Income", "(Average income attribute)", "HCPV", NA, NA, NA, NA, NA, NA)
)
diff5$level <- factor(
  diff5$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff5$estimate <- as.numeric(diff5$estimate)
diff5$std.error <- as.numeric(diff5$std.error)
diff5$z <- as.numeric(diff5$z)
diff5$p <- as.numeric(diff5$p)
diff5$lower <- as.numeric(diff5$lower)
diff5$upper <- as.numeric(diff5$upper)
diff5.plot <- ggplot(
  diff5, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff5.plot)

## MMs
res6 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ CPV_Group, level_order = "descending"
)
res6$level <- as.character(res6$level)
res6 <- rbind(
  res6, 
  c("HCPV", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCPV")
)
res6$level <- factor(
  res6$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res6$estimate <- as.numeric(res6$estimate)
res6$std.error <- as.numeric(res6$std.error)
res6$z <- as.numeric(res6$z)
res6$p <- as.numeric(res6$p)
res6$lower <- as.numeric(res6$lower)
res6$upper <- as.numeric(res6$upper)
res6.plot <- ggplot(
  res6, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = CPV_Group)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res6.plot)

diff6 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ CPV_Group, level_order = "descending"
)
diff6$level <- as.character(diff6$level)
diff6 <- rbind(
  diff6, 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCPV"), 
  c("HCPV - LCPV", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCPV")
)
diff6$level <- factor(
  diff6$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff6$estimate <- as.numeric(diff6$estimate)
diff6$std.error <- as.numeric(diff6$std.error)
diff6$z <- as.numeric(diff6$z)
diff6$p <- as.numeric(diff6$p)
diff6$lower <- as.numeric(diff6$lower)
diff6$upper <- as.numeric(diff6$upper)
diff6.plot <- ggplot(
  diff6, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff6.plot)


# COM
## AMCEs
res7 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ COM_Group, level_order = "descending"
)
res7$level <- as.character(res7$level)
res7 <- rbind(
  res7, 
  c("HCOM", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCOM")
)
res7$level <- factor(
  res7$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res7$estimate <- as.numeric(res7$estimate)
res7$std.error <- as.numeric(res7$std.error)
res7$z <- as.numeric(res7$z)
res7$p <- as.numeric(res7$p)
res7$lower <- as.numeric(res7$lower)
res7$upper <- as.numeric(res7$upper)
res7.plot <- ggplot(
  res7, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = COM_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res7.plot)

diff7 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ COM_Group, level_order = "descending"
)
diff7$level <- as.character(diff7$level)
diff7 <- rbind(
  diff7, 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Politics", "(Political attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Social", "(Social attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Population", "(Population attribute)", "HCOM", NA, NA, NA, NA, NA, NA), 
  c("HCOM - LCOM", "hlw_fsu", "amce", "Income", "(Average income attribute)", "HCOM", NA, NA, NA, NA, NA, NA)
)
diff7$level <- factor(
  diff7$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff7$estimate <- as.numeric(diff7$estimate)
diff7$std.error <- as.numeric(diff7$std.error)
diff7$z <- as.numeric(diff7$z)
diff7$p <- as.numeric(diff7$p)
diff7$lower <- as.numeric(diff7$lower)
diff7$upper <- as.numeric(diff7$upper)
diff7.plot <- ggplot(
  diff7, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff7.plot)

## MMs
res8 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ COM_Group, level_order = "descending"
)
res8$level <- as.character(res8$level)
res8 <- rbind(
  res8, 
  c("HCOM", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCOM")
)
res8$level <- factor(
  res8$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res8$estimate <- as.numeric(res8$estimate)
res8$std.error <- as.numeric(res8$std.error)
res8$z <- as.numeric(res8$z)
res8$p <- as.numeric(res8$p)
res8$lower <- as.numeric(res8$lower)
res8$upper <- as.numeric(res8$upper)
res8.plot <- ggplot(
  res8, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = COM_Group)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res8.plot)

diff8 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ COM_Group, level_order = "descending"
)
diff8$level <- as.character(diff8$level)
diff8 <- rbind(
  diff8, 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HCOM"), 
  c("HCOM - LCOM", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HCOM")
)
diff8$level <- factor(
  diff8$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff8$estimate <- as.numeric(diff8$estimate)
diff8$std.error <- as.numeric(diff8$std.error)
diff8$z <- as.numeric(diff8$z)
diff8$p <- as.numeric(diff8$p)
diff8$lower <- as.numeric(diff8$lower)
diff8$upper <- as.numeric(diff8$upper)
diff8.plot <- ggplot(
  diff8, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff8.plot)


# SS
## AMCEs
res9 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce", by = ~ SS_Group, level_order = "descending"
)
res9$level <- as.character(res9$level)
res9 <- rbind(
  res9, 
  c("HSS", "hlw_fsu", "amce", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "amce", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "amce", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "amce", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "amce", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HSS")
)
res9$level <- factor(
  res9$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res9$estimate <- as.numeric(res9$estimate)
res9$std.error <- as.numeric(res9$std.error)
res9$z <- as.numeric(res9$z)
res9$p <- as.numeric(res9$p)
res9$lower <- as.numeric(res9$lower)
res9$upper <- as.numeric(res9$upper)
res9.plot <- ggplot(
  res9, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = SS_Group)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(-0.4, 0.4) + 
  labs(y = "Estimated AMCE") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res9.plot)

diff9 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "amce_diff", by = ~ SS_Group, level_order = "descending"
)
diff9$level <- as.character(diff9$level)
diff9 <- rbind(
  diff9, 
  c("HSS - LSS", "hlw_fsu", "amce", "Distance", "(Distance attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_fsu", "amce", "Politics", "(Political attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_fsu", "amce", "Social", "(Social attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_fsu", "amce", "Administration", "(Administrative attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_fsu", "amce", "Population", "(Population attribute)", "HSS", NA, NA, NA, NA, NA, NA), 
  c("HSS - LSS", "hlw_fsu", "amce", "Income", "(Average income attribute)", "HSS", NA, NA, NA, NA, NA, NA)
)
diff9$level <- factor(
  diff9$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff9$estimate <- as.numeric(diff9$estimate)
diff9$std.error <- as.numeric(diff9$std.error)
diff9$z <- as.numeric(diff9$z)
diff9$p <- as.numeric(diff9$p)
diff9$lower <- as.numeric(diff9$lower)
diff9$upper <- as.numeric(diff9$upper)
diff9.plot <- ggplot(
  diff9, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.25, 0.25) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in AMCEs") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff9.plot)

## MMs
res10 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm", by = ~ SS_Group, level_order = "descending"
)
res10$level <- as.character(res10$level)
res10 <- rbind(
  res10, 
  c("HSS", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HSS")
)
res10$level <- factor(
  res10$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
res10$estimate <- as.numeric(res10$estimate)
res10$std.error <- as.numeric(res10$std.error)
res10$z <- as.numeric(res10$z)
res10$p <- as.numeric(res10$p)
res10$lower <- as.numeric(res10$lower)
res10$upper <- as.numeric(res10$upper)
res10.plot <- ggplot(
  res10, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
            ymax = upper, ymin = lower, shape = feature, color = SS_Group)) + 
  geom_hline(yintercept = 0.5, color = "gray60") + 
  geom_pointrange(position = position_dodge(width = 0.7), size = 0.2) + 
  scale_colour_manual(values = c("black", "grey60"), na.translate = FALSE) + 
  guides(shape = "none") + coord_flip() + ylim(0.25, 0.75) + 
  labs(y = "Estimated Marginal Mean") + 
  theme_minimal(base_size = 12) + 
  theme(plot.title = element_blank(), 
        axis.title.y = element_blank(), 
        legend.position = "bottom", 
        legend.title = element_blank())
plot(res10.plot)

diff10 <- cj(
  data2, 
  hlw_fsu ~ Distance + Politics + Social + Administration + Population + Income, 
  id = ~ no, estimate = "mm_diff", by = ~ SS_Group, level_order = "descending"
)
diff10$level <- as.character(diff10$level)
diff10 <- rbind(
  diff10, 
  c("HSS - LSS", "hlw_fsu", "mm", "Distance", "(Distance attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS - LSS", "hlw_fsu", "mm", "Politics", "(Political attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS - LSS", "hlw_fsu", "mm", "Social", "(Social attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS - LSS", "hlw_fsu", "mm", "Administration", "(Administrative attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS - LSS", "hlw_fsu", "mm", "Population", "(Population attribute)", NA, NA, NA, NA, NA, NA, "HSS"), 
  c("HSS - LSS", "hlw_fsu", "mm", "Income", "(Average income attribute)", NA, NA, NA, NA, NA, NA, "HSS")
)
diff10$level <- factor(
  diff10$level, 
  levels = c(
    "20 km or more", "5-20 km", "1-5 km", "less than 1 km", "(Distance attribute)", 
    "local economy", "public services", "job creation", "society", "(Political attribute)", 
    "70% agreed, 30% disagreed", "50% agreed, 50% disagreed", "30% agreed, 70% disagreed", "(Social attribute)", 
    "water charge", "public facility", "education", "elderly care", "childbirth and childcare", "none", "(Administrative attribute)",  
    "21,000", "20,000", "19,000", "(Population attribute)", 
    "4.2 million yen", "4 million yen", "3.8 million yen", "(Average income attribute)"
  )
)
diff10$estimate <- as.numeric(diff10$estimate)
diff10$std.error <- as.numeric(diff10$std.error)
diff10$z <- as.numeric(diff10$z)
diff10$p <- as.numeric(diff10$p)
diff10$lower <- as.numeric(diff10$lower)
diff10$upper <- as.numeric(diff10$upper)
diff10.plot <- ggplot(
  diff10, aes(x = reorder(level, as.numeric(feature) * -1), y = estimate, 
             ymax = upper, ymin = lower, shape = feature)) + 
  geom_hline(yintercept = 0, color = "gray60") + 
  geom_pointrange(size = 0.2) + coord_flip() + 
  ylim(-0.15, 0.15) + facet_wrap(. ~ BY, ncol = 1) + 
  labs(y = "Estimated Difference in Marginal Means") + 
  theme_minimal(base_size = 12) + 
  theme(axis.title.y = element_blank(), 
        legend.position = "none", 
        legend.title = element_blank())
plot(diff10.plot)
